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1. 


Summary  of  Achievements. 


Quantitative  prediction  of  the  structure,  properties,  and  chemical  reactivity  of  molecules 
by  computer  programs  using  first  principles  electronic  structure  methods  is  a  valuable 
complement  to  experimental  studies,  because  the  latter  are  often  difficult  or  impossible 
for  transient  reactive  species.  Unfortunately,  the  electronic  structure  methods  which 
achieve  quantitative  chemical  accuracy,  such  as  those  based  on  coupled  cluster  theory,  are 
currently  limited  in  applicability  to  small  molecules  of  less  than  about  ten  first  row 
atoms.  This  is  because  of  the  unphysical  increase  in  computational  cost  with  molecular 
size  in  standard  implementations.  In  this  Phase  I  SBIR  research,  Q-Chem  Inc.,  in 
collaboration  with  researchers  at  the  University  of  California,  Berkeley,  has  developed  a 
new-from-the-ground-up  approach  to  reformulating  these  proven  electron  correlation 
methods  in  terms  of  atom-centered  spatially  localized  quantities.  The  feasibility  and 
formal  advemtages  of  the  new  formulation  has  been  demonstrated  in  preliminary  work 
based  on  the  simplest  electron  correlation  method,  second  order  perturbation  theory. 

Trial  calculations  on  medium  sized  molecules  indicated  that  recovery  of  both  absolute 
electron  correlation  energy,  and  also  recovery  of  relative  energies  across  potential  surfaces 
is  excellent.  The  results  have  shown  the  potential  of  these  new  local  correlation  methods 
to  allow  the  study  of  much  larger  molecules  than  can  be  treated  conventionally,  with 
roughly  an  order  of  magnitude  improvement  being  possible. 
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2. 


Introduction. 


The  main  achievement  of  researchers  in  quantum  chemistry  over  the  past  three  decades 
has  been  the  development  of  a  rational  series  of  theories  of  molecular  wavefunctions 
which  form  successively  better  approximations  to  formally  exact  quantum  mechanics'. 
Each  approximation  may  be  viewed  as  a  “theoretical  model  chemistry”,  in  the  sense  that 
virtually  any  property  of  the  isolated  molecule  which  can  be  measured  experimentally  can 
be  calculated  within  the  chosen  framework.  The  goal  of  quantum  chemistry  is  the 
achievement  of  model  chemistries  which  are  sufficiently  accurate  to  permit  prediction  of 
molecular  structure,  properties  and  reactivity  in  cases  where  experiments  are  ambiguous 
or  unfeasible.  This  can  now  be  considered  virtually  accomplished  for  small  molecules"  '". 

As  a  result  of  this  systematic  development,  molecular  orbital  methods  have  been 
implemented  into  large  program  packages  which  have  been  eagerly  adopted  by  researchers 
in  industry  and  academia  alike.  There  are  now  roughly  ten  major  software  packages 
which  are  competing  for  sales  in  this  growing  market.  Applications  of  these  methods 
range  from  investigating  the  thermochemical  stability  of  chlorofluorocarbon  replacements 
to  calculations  of  the  partial  charges  on  small  biological  molecules  as  input  to  simulations 
for  computer-aided  drug  design.  Within  the  USA,  molecular  orbital  calculations  are  used 
as  a  tool  within  most  industrial  molecular  modeling  groups.  Part  of  maintaining  their 
international  competitiveness  is  to  ensure  that  such  groups  have  the  earliest  possible 
access  to  new  technical  advances,  as  has  historically  been  the  case  in  this  area. 

Yet  this  success  masks  an  absolutely  fundamental  problem,  which  must  be  addressed  if 
these  model  chemistries  are  ever  to  be  applied  to  molecules  containing  more  than  a  dozen 
or  so  atoms.  The  problem  is  that  the  computational  complexity  of  our  best  model 
chemistries  (for  example  based  on  the  Moller-Plesset  and  coupled  cluster  approaches) 
grows  in  an  unphysical  way  with  the  number  of  atoms.  While  electrons  interact  in 
pairwise  fashion,  the  complexity  of  the  coupled  cluster  approach  to  electron  correlation 
rises  as  the  sixth  power  of  the  size  of  the  molecule'''.  This  is  so  steep  a  dependence  that 
to  increase  the  size  of  molecule  by  a  factor  of  10  would  require  one  million  times  more 
computing  power  with  present  algorithms!  Even  computer  storage  requirements 
currently  rise  as  the  fourth  power  of  molecular  size,  so  ten  thousand  times  more  disk 
space  would  also  be  required. 

One  way  around  this  barrier  to  the  study  of  larger  systems  by  molecular  orbital  methods 
is  to  seek  entirely  new  theories.  For  this  reason,  much  attention  has  been  focused  on 
density  functional  theory  (DFT)  over  the  past  several  years.  Density  functional  methods 
scale  far  less  severely  with  molecular  size,  and  recover  the  proper  asymptotic  quadratic 
size  dependence  expected  from  pairwise  interactions''.  DFT  is  now  widely  employed  for 
quantum  mechanical  calculations  on  large  molecules.  However,  the  accuracy  of  DFT  is  all 
too  often  insufficient  for  chemical  purposes.  For  example,  barriers  to  chemical  reactions, 
which  are  critical  in  many  molecular  design  applications,  are  predicted  erratically  by 
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present  versions  of  DFT.  Hence  DFT  methods  are  not  proving  to  be  a  replacement  for 
accurate  Moller-Plesset  and  coupled  cluster  theory  calculations,  which  remain  the 
methods  of  choice  for  quantitative  calculations  when  they  are  computationally  feasible. 

The  other  alternative,  which  is  the  one  that  we  have  pursued  in  this  research,  is  to 
reformulate  high  level  molecular  orbital  theories  to  make  them  computationally  feasible 
for  studying  larger  molecules.  The  correlations  between  electrons  must  simplify  as  the 
separation  between  them  increases,  so  that  asymptotically  it  should  be  possible  to  obtain 
computational  complexity  that  increases  only  with  the  square  of  molecular  size, 
corresponding  to  pairwise  interactions.  The  length  scale  on  which  this  simplification 
occurs  can  be  estimated  very  roughly  from  the  length  scale  over  which  the  characteristic 
functional  group  properties  of  molecules  emerge.  This  is  on  the  order  of  only  a  few 
bondlengths,  which  suggests  that  for  molecules  containing  more  than  a  few  functional 
groups  there  will  be  tremendous  computational  savings  from  exploiting  the  natural 
distance  dependence  of  electron  correlations. 

We  have  developed  a  new-from-the  ground-up  tensor-based  formulation  of  electron 
correlation  during  the  Phase  I  period  which  opens  the  way  to  new  and  very  powerful 
“atoms-in-molecules”  (AIM)  methods  that  indeed  make  high  level  correlation  methods 
applicable  to  larger  molecules.  In  Sec.  3,  we  describe  the  underlying  basis  of  the  tensor 
formulation,  and  summarize  some  of  the  key  new  results  that  have  been  obtained. 

Perhaps  the  most  significant  one  is  the  fact  that  the  tensor  approach  permits  use  of 
nonorthogonal  atom-centered  functions  in  theories  of  electron  correlation.  In  Sec.  4,  we 
describe  two  proeedures  we  have  developed  for  obtaining  appropriate  sets  of  atom- 
centered  functions,  together  with  preliminary  results.  As  described  in  Sec.  5,  we  have 
combined  these  two  developments  to  obtain  computationally  efficient  AIM  treatments  of 
electron  correlation.  Any  standard  correlation  method  can  be  evaluated  by  our  AIM 
models.  We  have  demonstrated  their  accuracy  and  feasibility  via  a  pilot  implementation 
for  the  simplest  electron  correlation  method,  second  order  Moller-Plesset  theory. 

3.  A  Tensor  Approach  to  Electron  Correlation  Theory. 

In  the  Phase  I  period,  we  have  developed  a  straightforward  formalism  that  allows  theories 
of  electron  correlation  to  be  formulated  in  nonorthogonal  basis  sets'''.  Our  treatment 
makes  full  use  of  the  tensor  properties  of  the  operators  of  quantum  mechanics,  which 
permits  all  defining  equations  to  be  cast  in  a  manner  which  is  manifestly  invariant  to 
whether  or  not  the  basis  is  orthogonalized.  The  only  assumption  required  is  strong 
orthogonality  between  functions  spanning  the  occupied  subspace  and  functions  spanning 
the  virtual  subspace  of  the  one-particle  Hilbert  space.  Even  this  can  be  lifted  by  the 
introduction  of  projectors  into  the  equations.  The  tensor  approach  is  readily  employed 
using  algebraic  methods  in  either  first  or  second  quantization.  When  applied  to  standard 
correlation  methods,  such  as  coupled  cluster  theory  truncated  at  single  and  double 
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substitutions,  and  second  order  Moller-Plesset  perturbation  theory  we  obtain  equations 
which  are  invariant  to  larger  classes  of  transformations  than  existing  expressions. 


What  is  the  point  of  permitting  the  generalization  to  nonorthogonal  functions?  In  the  first 
instance  it  is  desirable  to  obtain  expressions  for  many-electron  theories  that  are  as  general 
as  possible.  One  might  imagine  that  given  such  equations,  new  and  useful  practical 
approximations  may  subsequently  arise.  For  example,  the  atomic  orbital  (AO)  single 
particle  basis  sets  ubiquitous  in  molecular  electronic  structure  theory  are  themselves 
nonorthogonal.  To  obtain  equations  that  are  directly  in  terms  of  this  given  basis  set,  there 
is  no  alternative  but  to  properly  build  nonorthogonality  into  the  formalism  from  the 
outset.  Since  these  functions  are  spatially  localized,  they  provide  a  natural  representation 
in  which  the  description  of  electron  correlation  via  many -body  theories  can  be  cast  in 
local  terms.  A  spatially  localized  treatment  of  electron  correlation''"’''"'  is  essential  for 
eliminating  the  unphysical  scaling  of  computational  cost  with  molecular  size  that  afflicts 
orthogonal  basis  formulations  of  many-body  theories. 

Why  was  a  formulation  of  this  kind  not  already  available?  All  standard  many  body 
methods  in  electronic  structure  theory,  be  they  first  quantization-based,  second 
quantization-based  or  diagrammatic,  implicitly  depend  upon  an  underlying  orthonormal 
single  particle  basis”'.  The  main  reason  for  employing  orthogonal  representations  is  that 
they  are  widely  believed  to  be  simpler.  In  particular,  the  usual  approach  to  obtaining 
equations  that  are  directly  in  a  nonorthogonal  basis  is  to  explicitly  carry  the  metric 
(overlap)  matrices  along  in  the  equations’'.  This  leads  to  equations  that  are  substantially 
more  cumbersome  in  the  nonorthogonal  basis  than  in  the  usual  orthogonal  representation. 
However,  with  the  full  use  of  tensor  algebra,  we  have  shown  it  is  possible  to  obtain 
expressions  that  are  general  to  nonorthogonal  basis  sets  with  relatively  little  additional 
complication. 


Space  considerations  prevent  us  from  giving  a  full  description  of  the  tensor  formalism 
here.  Instead,  we  shall  summarize  the  essential  concepts,  and  then  illustrate  the  form  of 
the  equations  that  arise  for  second  order  Moller-Plesset  theory.  The  equations  that  arise 
for  coupled  cluster  theory  are  also  straightforwardly  obtained,  and  will  be  explored  and 
implemented  as  part  of  our  Phase  II  research. 


Descriptions  of  tensor  analysis  are  available  in  many  textbooks’''”’'".  For  our  purposes  the 
essential  concepts  are: 

1)  Let  us  be  given  a  set  of  nonorthogonal  functions,  >}  The  way  in  which  these 

functions  transform  on  a  change  of  basis  is  termed  co variant,  and  covariant  indices  are 
denoted  with  subscripts.  The  overlap  matrix  of  these  functions  is  termed  the  metric: 

Sfiv  ~  ^fiv  ■ 
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2)  By  inverting  the  metric  matrix  (or  employing  a  generalized  inverse  if  necessary)  and 
acting  with  this  inverse,  =  (S~^)  ,  on  the  given  functions,  we  define  a  dual  basis, 

\  //XV 

V 

which  is  biorthogonal  to  the  original  functions:  ^0^  |^v  ^  •  Functions  in  the  dual 

basis  transform  in  an  inverse  way  to  functions  in  the  given,  covariant  basis,  and  indices 
with  this  property  are  termed  contravariant,  and  are  indexed  with  superscripts.  The 
action  of  the  metric  matrix  converts  indices  from  contravariant  to  covariant,  and  the 
inverse  metric  matrix  converts  from  covariant  to  contravariant. 

3)  By  defining  addition  and  contraction  of  tensors  appropriately,  tensor  algebra  yields 
results  which  are  by  construction  invariant  to  transformations  of  the  underlying  basis. 
Thus  contraction  of  a  pair  of  indices  (as  in  a  dot  product  of  vectors  to  yield  a  scalar)  is 
permitted  only  if  one  index  is  contravariant  and  the  other  is  covariant,  which  defines  the 
Einstein  summation  convention  (a  repeated  index,  one  co variant  the  other  contravariant, 
implies  summation).  Addition  of  tensors  is  only  possible  if  the  properties  of  their  indices 
match. 

4)  To  treat  the  electron  correlation  problem,  we  introduce  separate  sets  of 
nonorthogonal  functions  spanning  the  occupied  and  virtual  subspaces.  How  they  are 
obtained  is  discussed  in  Sec.  3.3.  Use  of  what  we  term  the  “natural  representation”  in 
which  many-electron  ket  quantities  contain  functions  in  the  covariant  basis,  and  many- 
electron  bra  vectors  contain  only  dual  basis  functions  leads  to  equations  for  electron 
correlation  methods  that  are  isomorphic  to  traditional  orthogonal  basis  derivations,  but 
substantially  more  general. 

5)  While  the  discussion  above  appears  to  involve  explicit  inversion  of  metric  matrices, 
which  appears  to  rule  out  the  use  of  linearly  dependent  representations,  this  is  not  in  fact 
the  case.  With  appropriate  use  of  the  theory  of  generalized  inverses’‘"',  and  the  singular 
value  decomposition,  it  is  quite  straightforward  to  employ  overcomplete  representations 
(as  we  in  fact  do  in  practice,  as  discussed  in  Secs.  4  and  5). 

We  give  an  example  of  the  tensor  formalism  applied  to  perhaps  the  simplest  useful 
treatment  of  electron  correlation.  This  is  the  many-body  perturbation  approach  first 
proposed  by  Moller  and  Plessef"'.  Revived  in  the  mid-1970’s,  this  approach,  carried  to 
second  order,  has  become  a  standard  tool  of  computational  quantum  chemistry'.  The 
usual  formulation  of  MP2  theory  involves  the  assumption  of  the  canonical  molecular 
orbitals  (so  that  the  zero  Hamiltonian  is  diagonal  in  the  many-electron  basis).  However, 
more  general  “noncanonical”  forms  of  MP2  theory  which  are  invariant  to  unitary 
transformations  within  the  occupied  and  virtual  spaces  were  introduced  in  the  context  of 
gradient  theory,  and  open  shell  problems’^''.  With  Pulay  and  Saebo’s  local  MP2  method''", 
a  yet  more  general  form  was  presented,  in  which  nonorthogonal  functions  were  employed 
to  span  the  virtual  function  space.  Our  tensor  formalism  yields  the  most  general  MP2 
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equations  yet  presented,  as  they  are  invariant  to  nonunitary  transformations  within  the 
occupied  and  virtual  spaces  individually. 

There  is  a  delightful  isomorphism  between  our  new  and  more  general  treatment  and  the 
conventional  orthogonal  basis  formulations.  The  first  order  equation  for  the  perturbed 
wavefunction  may  be  resolved  by  projection  with  singly  and  doubly  substituted  bra 
determinants  (in  the  dual  basis  of  course)  to  yield  equations  for  the  first  order  single  and 
double  substitution  amplitudes: 


(2a) 

(2b) 

These  equations  may  be  written  out  in  explicit  spin-orbital  notation  quite  easily,  using 
either  first  quantization  tensor-based  techniques  or  the  corresponding  second  quantization 
version.  Adopting  the  natural  representation  defined  above  for  all  tensor  quantities  we 
obtain: 

Fgl^-F/l^  +  Fi^^O  (3a) 

Fj-t"/  +  Fftff  -  -  F,hlj‘  +  nf  =  0  (3b) 

where  F  is  the  Fock  operator,  t  are  the  amplitudes  to  solve  for,  and  indices  are 
occupied,  a,b,c...  are  virtual,  and  of  course  subscripts  and  superscripts  have  their  tensorial 
meanings.  The  antisymmetrized  two  electron  integrals  are  explicitly  defined  in  the  natural 
representation  as: 

Finally,  the  second  order  expansion  of  the  energy  is  obtained  from  the  second  order 
expansion  of  the  eigenvalue  equation  followed  by  projection  with  the  bra  determinant: 

£«)  =(<i.|F<')|(r,")  +  r«))i.)=  fj;”  +1-0.“'’  (5) 

These  equations  constitute  the  most  general  form  of  the  MP2  equations  yet  presented: 
they  are  invariant  to  general  nonunitary  linear  transformations  within  the  occupied  and 
virtual  subspaces  respectively.  Within  the  framework  for  treating  linear  dependencies 
mentioned  above,  these  equations  can  even  be  solved  within  redundant  representations  if 
this  is  convenient.  These  equations,  despite  their  perturbational  nature,  must  be  solved 
iteratively  because  of  the  off-diagonal  Fock  terms  in  Eqs.  (3).  Finally,  Eqs.  (3)  can  be 
specialized  to  the  textbook  canonical  case  by  assuming  an  orthogonal  basis  (in  which  case 
the  covariant-contravariant  distinction  vanishes),  satisfying  the  Brillouin  theorem  (such 
that  first  order  singles  are  zero),  and  diagonalizing  the  occupied-occupied  and  virtual- 
virtual  blocks  of  the  Fock  matrix  (so  that  the  doubles  equations  uncouple). 

Via  the  use  of  metric  matrices,  it  is  straightforward  to  alter  the  tensor  character  of  the 
integrals  and  amplitudes  that  enter  the  above  equations.  For  example,  let  us  make  contact 
with  the  expressions  previously  presented  for  the  MP2  equations  in  a  representation 
where  the  virtual  subspace  was  treated  in  a  redundant  nonorthogonal  representation  while 
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the  occupied  subspace  was  described  in  a  noncanonical  but  orthogonal  basis''".  We  must 
make  essentially  only  one  change,  which  is  to  convert  to  the  fully  covariant  virtual 
representation  via  metric  matrices.  The  result  is  then  directly  comparable  with  the 
equations  that  form  the  starting  point  for  Pulay  and  Saebo’s  local  MP2  method''": 

=  0  (6a) 

(6b) 

In  their  local  MP2  method,  the  occupied  indices  are  restricted  to  an  orthogonal  basis,  so 
that  the  metric  matrices  for  the  occupied  space  are  simply  identity  matrices,  and  the 
distinction  between  covariant  and  contravariant  occupied  indices  vanishes.  They  also 
truncate  the  range  of  virtual  indices  associated  with  a  given  pair  of  occupied  indices  using 
a  spatial  proximity  condition''"’’^'''. 

To  summarize,  we  have  equations  for  second  order  Moller-Plesset  theory  (MP2)  and 
coupled  cluster  theory  which  are  more  general  than  any  previously  reported.  These 
equations  permit  the  use  of  local  nonorthogonal  functions  to  represent  both  occupied  and 
virtual  subspaces,  by  contrast  with  all  previous  formulations  which  have  required 
orthogonality  in  at  least  the  occupied  subspace.  Two  principal  issues  must  be  addressed 
in  order  to  begin  to  leverage  the  very  exciting  potential  of  this  new  formulation  for 
reducing  the  computational  complexity  of  electron  correlation  methods.  They  are: 

(1)  We  must  define  effective  techniques  for  obtaining  spatially  localized  sets  of 
functions  that  span  the  occupied  and  virtual  subspaces.  Our  solutions  to  this  problem  are 
discussed  in  Sec.  4.  This  bypasses  the  need  for  sets  of  localized  molecular  orbitals,  and 
the  problems  that  are  sometimes  associated  with  them. 

(2)  We  must  design  a  rational  hierarchy  of  local  correlation  models  that  can  be 
implemented  and  tested  with  the  functions  developed  above.  This  leads  to  our  “atoms- 
in-molecules”  (AIM)  methods  for  treating  electron  correlation,  described  in  Sec.  5,  which 
turns  out  to  have  unique  formal  and  practical  advantages  relative  to  previous  proposals 
for  local  MP2  methods. 

4.  The  Nonorthogonal  Orbitals. 

We  have  decided  to  focus  on  atom-centered  functions,  so  that  our  new  models  for  local 
correlation  can  be  based  on  atom-centered  selection  of  the  significant  pair  correlation 
amplitudes.  The  use  of  atom-based  functions  is  very  common  within  single  electron 
theories  of  electronic  structure  such  as  Hartree-Fock  theory  and  density  functional 
theory.  In  the  context  of  the  tensor-based  electron  correlation  formalism  developed  in 
Sec.  3,  we  require  sets  of  functions  spanning  the  occupied  and  virtual  subspaees.  The 
most  simple  approach  is  to  use  projected  atomic  orbitals:  projected  into  the  occupied 
space  and  the  virtual  spaee  respectively: 

=(f- 


(7a) 

(7b) 
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This  leads  to  functions  which  for  insulators  are  exponentially  localized,  provided  the 
system  under  study  has  a  band  gap.  This  condition  is  universally  satisfied  for  molecules. 
The  use  of  functions  of  this  type  is  satisfactory  for  the  virtual  space  and  has  been  used 
successfully  in  the  earlier  local  correlation  models  of  Pulay  and  Saebo''".  However 
projected  atomic  orbitals  are  several-fold  redundant  (and  hence  potentially  inefficient)  for 
the  occupied  space. 

We  have  hence  developed  alternative  treatments  of  the  occupied  space  which  reduce  its 
dimension  per  atom  to  that  of  a  minimum  basis  only  (instead  of  the  dimension  of  the  full 
basis  for  correlation  which  will  typically  be  3  to  10  times  larger).  It  has  been  long 
suspected  that  such  a  reduction  is  possible’'''"’’^''"'  with  little  loss  in  accuracy,  but  the  idea 
has  never  been  pursued  for  computational  advantage.  We  have  accomplished  the 
reduction  to  a  minimum  basis  in  two  alternative  ways,  which  each  have  different  special 
advantages.  We  shall  discuss  them  in  turn. 

Projected  Natural  Atomic  Orbitals  (PNAO*s).  This  is  an  exact  reduction,  based  on 
defining  natural  atomic  orbitals  from  atom-centered  blocks  of  the  one-particle  density 
matrix,  in  the  spirit  of  Weinhold’s  natural  atomic  orbitals’^'^.  Selecting  a  minimum  basis 
set  according  to  occupation  number  and  projecting  into  the  occupied  space  (via  the  analog 
of  Eq.  (7a))  yields  functions  that  exactly  span  that  space,  for  all  molecules  tested  (an 
extensive  series  of  compounds  comprised  of  first  row  atoms  and  hydrogen  atoms).  This 
set  of  localized  atom-centered  functions  gives  a  starting  point  for  correlation  treatments 
that  is  exactly  equivalent  to  the  conventional  choice  of  a  set  of  occupied  molecular 
orbitals.  The  extent  of  localization  is  directly  related  to  the  bandwidth  of  the  density 
matrix. 

Projected  Polarized  atomic  orbitals.  We  have  developed  an  alternative  reduction  in  which 
the  Hartree-Fock  calculation  that  defines  the  occupied  space  is  performed  in  a  molecule - 
optimized  atom-centered  minimum  basis  set.  We  term  these  functions  polarized  atomic 
orbitals  (PAO’s)  and  we  have  developed  a  pilot  program  that  determines  them  as  linear 
combinations  of  the  functions  included  in  a  much  larger  secondary  basis  set.  For  example, 
the  triple  zeta  6-3 1 IG**  basis  is  employed  to  define  the  molecule-optimized  minimal 
basis  of  PAO’s  in  the  preliminary  results  discussed  below.  Thus  given  a  (large) 

secondary  basis  the  (small)  basis  of  PAO’s  are  defined  via  an  atom- 

blocked  (let  PAO  index  a  belong  to  atom  A)  rectangular  transformation,  B: 

\<l>a)=  I  kKa 

HeA  ' 

where  the  transformation  coefficients  are  determined  as  part  of  the  variational  Hartree- 
Fock  calculation. 

The  PAO’s  appear  very  promising  as  a  method  for  performing  HF  (or  even  DFT) 
calculations  on  very  large  molecules  as  well  as  providing  atom-centered  functions  that 


(8) 
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span  the  occupied  space  (after  projection  with  the  density  operator).  While  the  Heirtree- 
Fock  energy  obtained  in  this  way  is  slightly  higher  than  that  in  the  full  (extended)  basis, 
this  difference  is  small  compared  to  normal  basis  set  incompleteness  effects.  More 
importantly,  the  quality  of  chemical  results  obtained  with  the  PAO  basis  is  vastly  higher 
than  that  of  a  conventional  minimal  basis,  and  approaches  that  of  the  full  secondary  basis. 

An  example  of  the  Phase  I  tests  of  P  AO’s  are  shown  in  Tables  1  and  2,  which  contains 
PAO  Hartree-Fock  and  PAO  second  order  Moller-Plesset  theory  calculations  on  a  small 
series  of  substituted  ethylenes,  probing  both  the  absolute  energy  of  the  most  stable 
conformer  (table  1),  and  the  relative  energies  of  the  stable  conformers  and  the  transition 
state  for  internal  rotation  (table  2).  As  these  molecules  are  all  conjugated,  they  represent  a 
stringent  test  for  the  performance  of  the  PAO  model  In  particular,  the  PAO  model  of 
bonding  approaches  exactness  in  the  limit  where  all  bonds  are  strictly  localized.  This  is 
clearly  not  the  case  in  these  conjugated  molecules,  where  the  size  of  the  rotational  barrier 
follows  directly  from  the  strength  of  the  conjugative  interactions  between  the  two 
unsaturated  groups. 


conformation 

HF/ 

HF/ 

MP2/ 

MP2/ 

PAO[6- 

311G**] 

6-31 IG** 

PAO[6- 

311G**] 

6-311G** 

Vi-CHO  (trans) 

-190.775497 

-190.812472 

-191.497914 

-191.504270 

Vi-OH  (trans) 

-152.904239 

-152.937832 

-153.473099 

-153.477798 

Vi-Vi  (trans) 

-154.929678 

-154.960122 

-155.592566 

-155.595937 

Vi-CFO  (trans) 

-289.656809 

-289.718733 

-290.600770 

-290.612309 

Table  1.  Hartree-Fock  and  second  order  Moller-Plesset  calculations  of  the 
absolute  energies  of  the  most  stable  conformer  of  a  series  of  conjugated 
ethylenes:  acrolein  (Vi-CHO),  vinyl  alcohol  (Vi-OH),  butadiene  (Vi-Vi),  and 
propenoyl  fluoride  (Vi-CFO).  All  calculations  are  at  HF/6-31G*  optimized 
geometries  and  the  absolute  energies  are  given  in  atomic  units  (Hartrees). 
Differences  between  using  the  polarized  atomic  orbitals,  and  full  use  of  the  6- 
31 IG**  basis  are  smaller  than  the  basis  set  truncation  errors.  Furthermore,  this 
difference  is  greatly  reduced,  by  a  factor  of  4  to  10,  at  the  MP2  level. 


conformation 

HF/ 

HF/ 

MP2/ 

MP2/ 

PAO[6- 

311G**] 

6-31 IG** 

PAO[6- 

311G**] 

6-311G** 

Vi-CHO  (trans) 

0.0 

0.0 

0.0 

0.0 
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Vi-CHO  (cis) 

2.162 

1.688 

1.610 

1.450 

Vi-CHO  (t.s.) 

6.649 

7.954 

7.691 

7.874 

Vi-OH  (syn) 

0.0 

0.0 

0.0 

0.0 

Vi-OH  (anti) 

2.777 

1.992 

2.131 

2.051 

Vi-OH  (t.s.) 

4.352 

4.574 

5.978 

5.573 

Vi-Vi  (trans) 

0.0 

0.0 

0.0 

0.0 

Vi-Vi  (gauche) 

3.199 

3.178 

2.730 

2.652 

Vi-Vi  (cis  t.s.) 

4.038 

4.098 

3.731 

3.699 

Vi-Vi  (perp.  t.s.) 

5.755 

6.058 

5.978 

5.770 

Vi-CFO  (trans) 

0.0 

0.0 

0.0 

0.0 

Vi-CFO  (cis) 

0.440 

0.198 

0.217 

0.093 

Vi-CFO  (t.s.) 

6.244 

7.604 

6.591 

6.674 

Table  2.  Hartree-Fock  and  second  order  Moller-Plesset  calculations  of  the 
relative  energies  of  the  stable  conformations,  and  transition  states  for  internal 
rotation,  of  the  same  series  of  substituted  ethylenes  given  in  table  1 .  All 
calculations  are  at  HF/6-31G*  optimized  geometries.  Relative  energies  are  given 
in  kcal/mol.  Differences  between  the  use  of  PAO’s  and  full  use  of  the  secondary 
basis  leads  to  changes  in  barrier  heights  which  are  generally  less  than  the 
estimated  uncertainty  of  +/- 1  kcal/mol  at  the  full  MP2/6-31 IG**  leveP’^. 

First  let  us  consider  the  results  in  table  1  briefly.  It  is  evident  that  the  difference  between 
PAO  and  full  energies  using  6-3 1 IG**  as  the  secondary  basis  is  only  the  order  of  several 
tens  of  milliHartrees.  By  comparison  the  total  energy  obtained  in  a  true  minimum  basis 
(STO-3G)  is  over  2.5  Hartrees  higher,  for  acrolein  as  a  typical  example!  The  PAO 
calculation  allows  us  to  quite  closely  approach  the  exact  energy  at  this  level  of  theory. 
Timiing  to  the  HF  values  of  relative  energies  shown  in  table  2,  we  see  that  the  difference 
between  PAO  relative  energies  at  the  HF  level,  and  results  obtained  with  the  full  6- 
3 1  IG**  basis  are  in  general  less  than  1  kcal/mol.  This  is  a  difference  which  is 
significantly  smaller  than  the  estimated  errors  in  the  calculation!  It  is  also  far  smaller  than 
the  energy  differences  between  the  PAO  calculations  and  the  full  6-3 1  IG**  calculations, 
which  is  on  the  order  of  20  kcaFmol.  Almost  entirely  all  of  this  difference  is  evidently 
common  to  the  various  isomers. 


Let  us  now  consider  the  similarities  and  differences  between  PAO  results  at  the  Hartree- 
Fock  level  of  theory,  and  PAO  results  at  the  correlated  MP2  level.  Inspection  of  the 
results  in  tables  1  and  2  suggests  the  following  main  points: 
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1)  From  table  1,  the  differenee  between  the  total  MP2  energy  in  the  secondary  basis, 
and  the  PAO-MP2  energy  is  lower  by  a  large  factor  (between  4  and  10)  than  the 
corresponding  Hartree-Fock  differences.  We  have  found  this  to  be  a  general  result,  which 
is  a  consequence  of  the  singles  contribution  to  the  MP2  energy  corresponding  to  an 
energy-weighted  steepest  descent  relaxation  of  the  density  matrix  obtained  in  the  PAO 
basis  towards  its  true  values. 

2)  As  a  result  of  this,  the  PAO-MP2  rotational  energy  barriers  approach  the  full  MP2 
rotational  barriers  somewhat  more  closely  than  the  PAO  Hartree-Fock  barriers 
approached  the  full  Hartree-Fock  barriers,  as  shown  in  table  2.  The  differences  between 
PAO  MP2  relative  energies  and  relative  energies  calculated  directly  in  the  secondary  basis 
is  less  than  the  estimated  uncertainty  of  about  +/-  1  kcal/mol  in  the  full  MP2 
calculations”^. 

3)  The  corresponding  differences  between  a  PAO-HF  reference  and  a  full  HF  reference 
in  an  even  better  correlation  treatment,  such  as  coupled  cluster  theory,  will  be  smaller  by 
at  least  another  equal  factor,  as  such  methods  can  easily  compensate  for  the  relatively 
small  difference  between  the  PAO  HF  reference,  and  the  full  HF  reference  in  the  basis 
used  for  the  correlation  treatment. 

The  use  of  PAO’s  offers  enormous  efficiency  advantages  for  the  treatment  of  electron 
correlation  at  the  second  order  Moller-Plesset  theory  level  This  advantage  derives  from 
the  fact  that  occupied  MO’s  can  be  directly  represented  in  terms  of  unprojected  PAO’s, 
which  constitute  only  a  minimal  basis.  Having  illustrated  that  the  PAO  model  is 
chemically  effective,  we  shall  develop  the  sophisticated  numerical  methods  necessary  to 
exploit  these  advantages  in  Phase  II.  The  availability  of  the  alternative  projected  atomic 
orbital  approach  ensures  that  we  can  also  offer  the  exact  solution  in  the  secondary  basis. 

5.  Atoms-in-Molecules  (AIM)  Methods  for  Electron  Correlation. 

The  developments  described  in  the  previous  two  sections  constitute  a  formally  exact  re¬ 
expression  of  standard  theories  of  electron  correlation.  In  the  traditional  formulation,  the 
nonredundant  canonical  molecular  orbitals  are  employed  to  express  many-electron 
theories.  The  canonical  MO’s  are  in  general  fully  delocalized  throughout  the  system.  In 
our  re-expression  of  electron  correlation  theories,  redundant  sets  of  atom-centered 
funetions  are  used  to  develop  the  working  equations.  The  question  addressed  in  this 
section  is  how  we  exploit  the  spatial  localization  of  these  atom-centered  functions  to 
develop  new  local  correlation  models  that  involve  drastically  reduced  computational 
complexity.  The  validity  of  strong  local  truncations  was  first  demonstrated  by  Pulay’'’^'. 
Our  tensor-based  formalism  will  permit  comparable  models  to  be  formulated  directly  in 
terms  of  atom-centered  truncations,  with  signifieant  formal  and  practical  advantages. 


12 


First  we  discuss  a  completely  general  procedure  for  solving  the  MP2  equations  in  a 
subspace  of  the  full  four-index  set  of  pair  correlations.  Let  us  initially  make  the  two- 
electron  integrals  in  the  fully  covariant  representation, 

(1)  Definition  of  the  truncated four-index  ket  space.  We  discard  a  large  number  of  the 
fully  covariant  integrals  based  on  the  models  described  below,  leaving  us  with  a  reduced 
set  in  which  the  indices  a,b,i,j  are  no  longer  independent.  Instead  we  can  view  the 
retained  abij  as  defining  a  single  linear  index,  (abij).  In  tensor  terms,  this  truncation  has 
reduced  the  rank  of  the  four  index  space  from  4  to  1 .  Effectively  we  have  selected  a 
subset  of  the  ^<^(ijab^  within  which  we  shall  satisfy  the  equations.  We  will  surround 

the  compound  index  with  brackets  to  indicate  that  the  individual  indices  are  no  longer 
independent. 

(2)  The  definition  of  the  metric  of  the  truncated  four-index  space.  Have  admitted  that  in 
the  truncated  space,  we  can  no  longer  permit  the  indices  a,b,i,j  to  be  independent,  we  can 
no  longer  define  the  dual  4  index  space  in  terms  of  simple  products  of  the  inverses  of  the 
single  index  metrics.  For  example,  simply  defining  the  dual  space  of  truncated  bra 
functions  in  that  manner  will  lead  to  local  correlation  models  in  which  the  correlation 
equations  are  not  solved  in  a  well-defined  subspace.  Instead  we  first  properly  define  the 
metric  associated  with  our  truncated  ket  space.  This  takes  the  general  form: 

^(ijab)(klcd)  ~  SikSjlSacSbd 

If  we  need  to  transform  a  truncated  tensor  of  character  into  character  then 
the  metric,  (9),  must  be  inverted  either  implicitly  or  explicitly: 

^(klcd){ijaS).  ikied) 

(3)  Solution  of  the  truncated  equations.  We  wish  to  solve  the  tensor  MP2  equations 
within  the  subspace  selected  in  Step  (1)  above,  with  metric  characterized  in  Step  (2). 

This  leads  us  to  the  following  explicit  form  of  the  truncated  equations,  in  terms  of  the 
corresponding  set  of  truncated  amplitudes: 

^^acEbd"^  Sac^bd}^ikSjl  ^  ~  ~  ^(ijab) 

This  set  of  equations  explicitly  defines  the  form  of  the  coefficient  matrix  in  the  truncated 
4-index  space.  The  equations  can  be  solved  by  standard  conjugate  gradient  or  DIIS’^’^" 
techniques,  although  we  note  that  as  written  here,  the  form  of  the  residual  on  each 
iteration  is  not  tensorially  consistent  with  the  form  of  the  amplitudes  that  we  must 
update.  Conversion  is  accomplished  via  solution  of  Eqs.  (10)  (which  need  not  be  exact 
for  this  purpose).  Finally,  given  the  converged  amplitudes  in  the  truncated  space,  we 
obtain  the  MP2  energy  in  the  form: 

(12) 

When  the  correlation  equations  are  solved  in  a  truncated  4-index  space  in  the  above 
manner,  we  aehieve  solution  of  the  equations  within  a  well-defined  subspace  of  the  full 
set  of  double  excitations.  By  contrast,  if  we  did  not  employ  the  proper  metric  associated 


(9) 

(10) 

(11) 
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with  the  truncated  rank  1,  4  index  space,  then  the  solution  lacks  this  important  model 
property.  Hence  our  new  methods  are  rigorously  invariant  to  linear  transformations  that 
mix  the  retained  excitations. 

Within  the  Phase  I  project,  we  have  formulated  a  series  of  truncations  that  are 
increasingly  accurate.  Since  the  truncations  are  based  on  our  atom-centered  functions 
spanning  occupied  and  virtual  spaces,  we  refer  to  these  as  “Atoms  in  Molecules”  (AIM) 
models  of  localized  electron  correlation  effects.  We  restrict  the  admissible  correlation 
amplitudes,  ,  via  the  following  models.  We  note  that  models  such  as  these  are 
entirely  new,  because  they  could  not  be  constructed  without  the  new  developments 
described  in  Secs.  3  and  4. 

(1)  Covalent  two-center  truncatiorf^"’.  Compound  indices  (abij)  are  retained  if: 

•  indices  i  and  a  belong  to  the  same  center,  and  j  and  b  belong  to  the  same  center.  These 
interactions  correspond  to  direct  dispersion  interactions  at  long  range,  and  decay  as 
inverse  cube  of  length  with  the  separation  of  the  sites  (assuming  local  neutrality). 

•  indices  i  and  b  belong  to  the  same  center,  and  j  and  a  belong  to  the  same  center.  These 
are  the  exchange  interactions  that  complement  the  direct  correlation  effects  above, 
ensuring  a  balanced  approximation.  These  interactions  (in  insulators)  decay  exponentially 
with  the  separation  of  the  separation  of  the  sites,  and  thus  for  large  molecules  (in  the 
range  of  tens  of  atoms)  it  becomes  profitable  to  apply  cutoffs  to  these  terms.  This  model 
is  in  a  sense  the  AIM  analog  of  the  original  Pulay  and  Saebo  proposal  for  a  local  ansatz''". 

(2)  Full  two-center  model.  The  most  sophisticated  two-center  model  is  to  retain 
all  compound  indices  where  no  more  than  two  centers  are  involved  in  all  four  indices. 

Only  the  dispersion-like  interactions  (identified  in  model  (1))  are  truly  long-range,  and 
thus  cutoffs  can  be  profitably  applied  to  the  model  approximation.  This  model  is  exact 
for  correlation  of  bond  orbitals  that  are  strongly  localized. 

The  covalent  two-center  model  can  be  shown  to  have  several  very  beautiful  formal 
properties  which  make  it  particularly  attractive  as  a  target  theory.  In  particular,  we  have 
shown  that: 

(1)  The  metric  is  a  separable  product  (because  of  its  diatomic  nature).  This  means  that 
iterating  to  invert  the  full  metric  is  not  necessary.  Instead  we  have  explicitly  inverted  the 
mono-atomic  metric  via  a  singular  value  decomposition,  and  then  formed  the  inverse  of 
the  diatomic  metric  as  a  tensor  product.  This  entirely  eliminates  the  necessity  to  solve 
Eq.  (10),  which  is  the  first  iterative  step  in  the  general  local  scheme  described  above. 

(2)  It  is  possible  to  recanonicalize  the  MP2  equations,  (1 1),  by  exploiting  the  product 
nature  of  the  metric.  In  particular  the  coefficient  matrix  that  defines  Eq.  (11)  separates 
into  the  sum  of  two  products  of  lower  rank  matrices,  one  of  which  is  the  fundamental 
mono-atomic  metric  matrix,  and  the  other  of  which  is  a  mono-atomic  excitation-like 
matrix.  The  former  we  have  already  diagonalized  in  the  first  step  above.  The  latter  can  be 
diagonalized  in  the  basis  of  eigenvectors  of  the  metric,  which  makes  the  full  MP2 
coefficient  matrix  diagonal.  As  a  result,  the  AIM-MP2  method  is  entirely  noniterative,  just 


14 


like  conventional  MP2,  and  entirely  unlike  previous  local  MP2  methods.  The 
computational  requirements  for  these  two  diagonalization  steps  are  strictly  cubic  with  the 
size  of  the  molecule. 

As  is  evident  from  the  description  above,  there  are  tremendous  practical  advantages  to 
employing  the  first  AIM-MP2  model,  as  it  can  be  made  entirely  noniterative.  Therefore 
we  have  carefully  investigated  its  behavior  with  respect  to  recovery  of  the  exact  MP2 
correlation  energy  to  see  if  it  is  a  suitable  model  for  practical  chemical  calculations.  For 
the  preliminary  investigations  completed  in  this  Phase  1  work,  we  employed  projected 
atomic  orbitals  for  both  occupied  and  virtual  subspaces  to  permit  direct  comparison  with 
conventional  canonical  MP2  calculations.  A  pilot  code  was  developed  which  was  not  in 
any  way  optimized,  hence  we  cannot  report  elapsed  times  at  this  stage.  Rather  the 
purpose  of  these  calculations  is  to  assess  the  quality  of  the  chemical  information  that  the 
new  AIM  models  are  capable  of  yielding.  We  know  from  straightforward  numerical 
analysis  of  the  new  models  that  they  are  capable  of  being  applied  to  very  large  systems, 
but  achieving  that  goal  is  only  worthwhile  once  it  has  been  established  that  the  results  are 
of  sufficient  quality  to  be  useful. 

Tables  3  and  4  below  summarize  the  results  of  a  series  of  test  calculations  we  have 
performed  on  conjugated  molecules  similar  to  those  used  to  test  the  polarized  atomic 
orbital  model  in  Sec.  4.  Table  3  reports  the  recovery  of  the  absolute  MP2  correlation 
energy,  and  it  is  evident  that  AIM-MP2  is  capable  of  yielding  roughly  99.9%  of  the 
correlation  energy  for  these  systems.  This  is  a  very  high  fractional  recovery  of  the 
electron  correlation  energy,  which  is  sufficiently  good  that  we  should  expect  the  recovery 
of  relative  correlation  energies  to  be  excellent  also.  We  note  parenthetically  that  the 
problem  of  relative  correlation  energies  is  much  tougher  than  absolute  correlation  energies 
because  we  are  concerned  with  small  energy  differences  (perhaps  on  the  order  of  kcal/mol) 
between  quite  large  correlation  energies  (on  the  order  of  Hartrees,  or  between  100  and 
1 000  kcal/mol).  To  obtain  less  than  10%  error  in  a  relative  energy  of  that  magnitude 
requires  a  differential  error  of  less  than  about  0.03%  between  the  points  of  interest  on  the 
potential  surface,  which  is  challenging  indeed. 


molecule 

E(HF)/a.u. 

dE(MP2)/a.u. 

dE(AIM-MP2) 

%  correlation 

Vi-C2H3 

-154.960121 

-0.560776 

-0.560541 

99.96 

Vi-CFO 

-289.718733 

-0.799584 

-0.799553 

99.99 

Vi-CHO 

-190.812152 

-0.604500 

-0.604461 

99.99 

Vi-NH2 

-133.103964 

-0.466910 

-0.466784 

99.97 

Vi-N02 

-281.578165 

-0.853139 

-0.853124 

99.99 

Vi-OH 

-152.937831 

-0.483513 

-0.483439 

99.98 

15 


Ph-CHO 

-343.511145 

-1.140552 

-1.138509 

99.82 

Ph_NH2 

-85.800480 

-1.002452 

-1.000878 

99.84 

Ph-OH 

-305.634634 

-1.017709 

-1.016479 

99.88 

Table  3.  Recovery  of  the  canonical  MP2  correlation  energy  by  the  local  AIM- 
MP2  method.  Calculations  are  using  the  6-31 IG**  basis  set,  at  the  lowest 
energy  conformation  of  each  molecule,  which  was  optimized  at  the  HF/6-3 IG* 
level.  The  final  column  is  the  percentage  of  exact  correlation  energy  recovered 
by  the  AIM-MP2  model,  which  is  typically  in  the  vicinity  of  99.9%. 


molecule 

A(HF)/kcal  mof' 

A(MP2)/keal  mof’ 

A(AIM-MP2) 

Vi-C2H3 

6.058061 

0.313 

0.270 

Vi-CFO 

7.604045 

0.955 

0.960 

Vi-CHO 

8.006900 

0.106 

0.109 

Vi-NH2 

6.477682 

0.745 

0.789 

Vi-N02 

6.955837 

1.834 

1.787 

Vi-OH 

4.573220 

0.975 

0.958 

Ph-CHO 

8.809472 

0.527 

0.722 

Ph-NH2 

3.746802 

0.510 

0.665 

Ph-OH 

2.570867 

0.788 

0.723 

Table  4.  Differential  recovery  of  the  canonical  MP2  correlation  energy 
contribution  to  rotational  energy  barriers  by  the  local  AIM-MP2  method. 
Calculations  are  performed  using  the  6-31  IG**  basis  set,  at  geometries  which 
were  optimized  at  the  HF/6-3  IG*  level.  All  energies  are  reported  in  units  of 
kcal/mol.  It  is  evident  that  the  AIM-MP2  method  recovers  these  very  small 
energy  differences  with  errors  of  less  than  10%.  As  discussed  in  the  text,  it  is 
not  clear  that  these  errors  are  to  be  regarded  as  deficiencies  relative  to 
experiment. 


Table  4  contains  the  results  of  relative  correlation  energy  calculations  for  the  rotational 
barriers  in  the  conjugated  molecules.  It  is  clear  that  AIM-MP2  recovers  the  exact  results 
to  quite  aeeeptable  accuracy.  On  formal  grounds,  we  expect  these  deviations  to  generally 
decrease  as  the  basis  set  is  further  improved.  Furthermore,  on  formal  grounds  it  is  not 
obvious  that  the  “exact”  results  are  to  be  preferred.  As  has  been  illustrated  previously  by 
Pulay  and  co-workers  for  the  water  dimer,  and  Friesner  and  co-workers  for 
conformational  energy  differences  in  their  local  method,  basis  set  superposition  error 
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(BSSE)  affects  the  exact  results  far  more.  The  principal  conclusion  we  draw  from  tables  3 
and  4  is  that  the  quality  of  our  simplest  AIM  model  is  sufficiently  high  that  we  are  able 
to  exploit  its  desirable  properties,  and  we  are  not  forced  at  present  to  pursue  the  most 
general  4  center  local  model.  This  is  pleasing  both  on  physical  grounds  (for  the  model  is 
simpler),  and  on  practical  grounds  (the  model  is  much  faster  and  more  efficient). 

6.  Epilog. 

With  the  new  Atoms-in-Molecules  (AIM)  approach  for  reducing  the  complexity  of 
electron  correlation  methods  in  hand,  the  stage  is  set  for  the  Phase  II  effort.  This  is  the 
extension  and  efficient  implementation  of  AIM  methods  for  the  highest  level  electron 
correlation  treatments  in  routine  use  for  small  molecules:  coupled  cluster  theory.  The 
result  will  be  a  unique  computational  tool  that  permits  routine,  highly  accurate  treatment 
of  problems  that  are  in  the  range  of  an  order  of  magnitude  larger  than  is  possible  at 
present,  using  methods  whose  accuracy  is  already  accepted  and  established  beyond 
doubt! 
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